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CONTRACTOR REPORT 


A COMPUTER CODE FOR THREE-DIMENSIONAL INCOMPRESSIBLE FLOWS 
USING NONORTHOGONAL BODY-FITTED COORDINATE SYSTEMS 

INTRODUCTION 


With the currently increasing computer capability and various flow solvers 
developed, numerical simulations of three-dimensional incompressible flow problems 
using Reynolds -average Navier-Stokes equations are now becoming more feasible in 
many engineering design and analysis applications. In many real world flow problems, 
the boundary geometries are complex such that it is more accurate to describe the 
geometries using body-fitted coordinate (BFC) systems. Especially for internal flow 
problems with complex geometries such as those of the hot gas manifold (HGM) of the 
Space Shuttle Main Engine (SSME), the use of nonorthogonal BFC systems for numeri- 
cal solutions can be beneficial in many aspects. It is not only the boundary geomet- 
ries that can be represented more closely using BFC systems, but also grid-refined 
solutions can be obtained without increasing an excessive amount in computer memory. 
In addition, once a particular flow problem has been set up, the redesign or opti- 
mization process of the boundary shapes can be performed very easily using BFC 
systems. 

Several numerical methods [1,2, 3, 4, 5, 6] has been developed for solving the 
incompressible Navier-Stokes equations in 3-D BFC systems. The main difference 
between these methods lies in the way of finding a pressure field such that the 
flowfield can be as close to divergence -free as possible (i.e. to satisfy the mass con- 
servation equation) . This is the main feature and difficulty of solving the incom- 
pressible flow problems. Numerical methods of References 1, 2 and 3, for instance, 
have employed the pseudocompressibility approach and time-iterative scheme to gen- 
erate the pressure field so that the continuity equation is satisfied when a steady 
state solution is reached. In these methods, artificial smoothing techniques must be 
used to obtain a strong coupling between the velocity and pressure fields. Methods 
of References 4, 5 and 6, on the other hand, have utilized a successive pressure- 
velocity correction scheme by using a Poisson’s equation for pressure correction 
derived approximately form the continuity and momentum equations. For these latter 
methods, grid staggering between the velocity vectors and the pressure nodes must 
be used to ensure stability of the numerical solutions. 

There are several possible methods of grid staggering associated with different 
features in solving the pressure correction equation. These grid staggering methods 
were discussed in Reference 6, from which one of the methods was shown to be the 
most promising arrangement (i.e. with the velocity vectors located at the faces of a 
volume which contains the pressure and other scalars at its center). But, this 
method has one drawback, that the velocity components are solved using different 
control volumes. It is for this reason that a grid staggering system similar to the 
one used by Vanka et al. [4] is developed in the present study. The present 
method of grid staggering and pressure correction equation that was described by 
Vanka [4[ and Maliska [6] . Also, using the present method, the same control volume 
is used for the velocity components and scalar quantities. 



In the following sections, basic elements for establishing the present computer 
code for solving the curvilinear Navier- Stokes equations in three-dimensional space 
(CNS3D) will be described. These are followed by a series of standard numerical 
examples used to evaluate the accuracy and efficiency of the present numerical 
method. The numerical examples include laminar flow driven-cavity problem, cases 
of laminar and turbulent flows over backward- facing steps, and 3-D laminar flows 
inside a 90-deg-bend square duct. Applications of the present code to the internal 
flow problems of SSME will be included in future publications. 

A user’s guide to the present CNS3D code is provided in Appendix A. 
Appendix B contains a list and definitions of all the major fortran symbols used in 
the computer program which is listed in Appendix C . 


TRANSFORMATION OF THE EQUATIONS OF MOTION 


For incompressible Newtonian fluid, the continuity, momentum and energy 
equations can be written as: 


U t + E x +F y +G z = S <» 


where (x,y,z) represent the Cartesian coordinates, and 
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e = the internal energy per unit mass = C y T for perfect gas 


Q = energy added per unit volume 
k T = thermal conductivity of the fluid 


Equation (1) is transformed to a general curvilinear coordinate system (£, n, 
O, which results in equation (2). 
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The transformation coefficients, £ x > £ , ? z , n x> n y , n z , C x , Cy, and c z > are com- 
puted numerically using second order central differencing. In the transformed domain, 
the grid sizes (i.e., A £ , An, and Ac) are set to be unity. This simplifies the cal- 
culation of the transformation coefficients. 

For turbulent flow computations, the present code has employed the standard 
k-e turbulence model [7] to provide the turbulent eddy viscosity The standard 

k-e turbulence model (which consists of a turbulent kinetic energy equation, k- 
equation, and a turbulent kinetic energy dissipation rate equation, e-equation) is 
given as: 


(Pk) t + 


(p e) t + 
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where the effective viscosity is calculated from: 
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and the turbulent kinetic energy production term, P f , is defined as: 
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The turbulence model constants are: 


C - 0.09 

P 

C x = 1.44 


a. = 1 . 0 , a-1.3 

C 2 - 1.92 . 
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Also, the molecular viscosity u in equation (1) is replaced by the effective viscosity 
y eff ^ or turbulent flow cases. 

In order to save the computational efforts, the widely used wall function 
approach [ 8] is employed to provide the near wall boundary conditions for the momen- 
tum and energy equations and the k-e turbulence model. This approach avoids the 
requirement of integrating the governing equations up to the wall which requires a 
large number of additional grid points near the wall. 

Equations (2), (3), and (4) form a closed set of nonlinear partial differential 
equations governing the fluid motion. This set of equations are to be solved by 
means of finite difference approximations which are performed in the transformed 
domain. For treating the convection terms, the hybrid scheme [9] is employed for 
simplicity (although other more elaborate schemes such as central differencing plus 
artificial dissipation scheme, QUICK scheme, or skew upwind differencing scheme, etc. 
can be implemented [10]). These are described in the following sections. 


DISCRETIZATION OF THE EQUATIONS OF MOTION 


In this section, finite difference approximations are used to discretize the 
governing equations, equations (2), (3), and (4). Second-order central differencing 
is used for the diffusion terms and the source terms. The hybrid differencing scheme 
[9] is employed to approximate the convection terms in the governing equations. The 
finite difference discretizations are performed in the transformed domain. The solu- 
tion procedure for the discretized equations using a velocity-pressure correction 
algorithm (SIMPLE-C) of References 11 and 12 will be described in the next section. 


The governing equations of motion can be represented by the following model 
transport equation in which <f> denotes all the dependent variables respectively and r 
is the diffusion coefficient.. 
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Discretization of equation (5) is performed using finite difference approximations in 
the transformed domain. The second order central differencing is used for approxi- 
mating the diffusion terms. For the convection terms, the hybrid differencing scheme 
[9] is employed (i.e., using central differencing for cell Peclet number less than or 
equal to 2 and switching to upwind differencing when the cell Peclet number is 
greater than 2) . The finite difference equation is arranged by collecting terms 
according to the grid nodes around a control volume as shown in Figure 1. The final 
expression is given by equation (6) in which A represents the link coefficients 
between grid nodes P, E, W, N, S, T, B, NE, NW, NT, NB, SE, SW, ST, SB, ET , 
EB , WT, and WB as shown in Figure 1. 
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( 6 ) 


The subscript ° denotes the solution at the previous time level. A fully implicit 
formulation is employed for solving the time dependent transient problems. 



Figure 1. Three-dimensional grid structure and labeling around a grid node P. 


6 


Thus, the nonlinear equations of motion are approximated by a system of linear 
algebraic equations which have the form of equation (6). Only one program subrou- 
tine is designed to calculate the link coefficients and the source terms. The number 
of algebraic equations depends on the number of interior grid points. For a grid 
size of 10 x 10 x 10 the number of algebraic equations to be solved would be around 
512. This large system of equations are preferred to be solved by some iterative 
methods, such as Gauss-Seidel iteration, line- underrelaxation method [13] or Stone's 
method [14], etc., rather than using direct methods such as Gaussian elimination 
method. Only a few (6 to 10) iterations through the whole computational domain are 
needed and a complete convergence of the system of algebraic equations is not 
required. Since equation (6) is only a linearized version of the governing equations 
which are nonlinear and coupled in nature, solutions of the equations of motion must 
be obtained through global iterations among the equations. A tentative solution to 
equation (6) will not affect the final results significantly. On the other hand, if too 
many iterations are used to get a better solution of equation (6), then a great deal 
of computing time would be virtually wasted. However, the above argument can not 
be applied when the pressure correction equation (which will be derived in the next 
section) is solved. Since during each global iteration it is desirable to retain a 
divergence -free velocity field, better solution of the pressure correction equation 
would in effect promote the convergence of the whole numerical scheme. Therefore, 
more iterations are usually used to solve the pressure correction equation. 


SOLUTION PROCEDURES 


The governing equations used in the present analysis are nonlinear and strongly 
coupled. Iterative procedures are employed to drive the equations to a converged 
solution. It is particularly important for incompressible flow to make the flow field 
satisfy the continuity equation and the momentum equations at the same time. This 
requires a correct pressure field associated with a divergence- free velocity field. 

A velocity-pressure correction procedure is developed in the present study to drive 
the pressure field and the velocity field to be divergence free. This kind of pro- 
cedure requires grid staggering between the velocity components and the locations 
where the pressure is estimated and stored such that the velocity field and the 
pressure field will not be uncoupled. 

In the present study, staggering grid systems as shown in Figure 2 (for 2-D 
case) are used. The velocity components, u and v, are solved and stored at the 
grid nodes and the pressure, p, is located at the corners of the control volume of 
u and v. In this way, solutions of u and v can be solved using the same control 
volume and coupling between u, v and p can also be enforced. To estimate the 
pressure field, a pressure correction equation is derived approximately from the dis- 
cretized momentum and continuity equations. The velocity and pressure fields are 
then corrected using the solutions of the pressure correction equation. 

First, the finite difference momentum equations (for u, v and w) can be written 
as: 
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u T 
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u, v AND SCALARS 
PRESSURE 


Figure 2. Locations where the variables are stored 
(staggering grids are used). 
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P 
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w 



(8b) 


(8c) 


where u*, v*, w*, and p* represent the solutions of equations (8a) and (8b). To 
satisfy the continuity equation the velocities and pressure are corrected according to 


the following relations: 


U = U* + U T 

(9a) 

v — V* + V T 

(9b) 

w = w* + w 1 

(9c) 

p = p* + pt 

(9d) 
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A new set of momentum equations can be constructed approximately using the 
divergence- free flow field, u, v, w, and p: 
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w 


(10a) 
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(10c) 


By subtracting equations (8a) through (8c) from equations (10a) through (10c), 
respectively, the following equations result: 
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According to SIMPLE-C algorithm [11], equations (11a) through (11c) are rearranged 
to be: 


<A p U - £ A ( u ) u p - = £ Aj u (Uj' - u p ') - P x ' 
i i 
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i i 
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The first terms on the right-hand side of equations (12a) through (12c) are neglected 
to simplify the formulation. Thus, 






i 


(13a) 


(13b) 


(13c) 


Using the decompositions of equations (9a) through (9c), the continuity equation can 
be written as: 


u + v + 

x y 


w 


= (u * + V * + W *) + (u ' + V, 


w z’) = 0 


(14) 


Substituting equations (13a) through (13c) into equation (14), the following pressure 
correction equation can be obtained: 


-[(D P ') + (D P ’) + (D P ') ] = -(u * + v * + w *) 

v u x y x v v y 'y v w z z J v x y z 


(15) 


Equation (15) is a Poisson’s equation with the source term equal to the local diver- 
gence of the flow field. To enforce the coupling between the velocity and pressure 
fields, the source term of equation (15) is first evaluated at the control volumes 
centered between the velocity nodes as shown in Figure 3. An averaged source term 
is then calculated at the cell center of p node for solving equation (15). In this way, 
the difficulties in solving the pressure correction equation, as described by Vanka [4] 
and Maliska [6], are eliminated. Coupling between the velocity and the pressure field 
is also assured. 

According to the above analyses, the present numerical method contains the 
following solution steps: 

1) Guess initial velocity and pressure field. 

2) Solve for the velocity field using equations (8a) through (8c) . 

3) Solve for other scalar transport equations. 

4) Solve the pressure correction equation, equation (15). 
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Figure 3. Control volumes where the mass conservation is evaluated 
for solving the pressure correction equation. 

5) Correct the velocity and pressure fields using equations (13a) through 
(13c) and equation (9d). 

6) Go back to step (2) until solution converges. 

A converged solution is obtained when the following criterion is met: 


Error = ( Au 


max 


+ A v 


1 max 


+ Aw ) / U o + P' /pU f 
max ref max ref 


i 3 x 10 


-4 


where Au, Av, and Aw represent velocity changes during each iteration due to the 
solutions of the momentum equations. 

In solving the momentum equations in step (b) above, underrelaxation factor 
of about 0.6 is recommended. With this, A p 's in equations (8a) through (8c) are 

modified according to the underrelaxation factor. For the correction of velocity field, 
no underrelaxation is required. But the correction of pressure field should be under- 
relaxed slightly (around 0.9) when the grid nonorthogonality is strong. This is 
different from that suggested by References 11 and 12 (which recommend no under- 
relaxation for pressure correction) . 
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NUMERICAL EXAMPLES 


In this section, several numerical examples are employed to demonstrate the 
efficiency and accuracy of the present numerical method. To serve this purpose, 2-D 
and 3-D, laminar and turbulent flow cases are included. These cases are: (a) 2-D 

laminar driven square-cavity flows; (b) 2-D laminar flows over a backward- facing 
step; (c) 2-D turbulent flows over a backward- facing step; (d) 3-D developing 
laminar flow inside a 90-deg-bend square duct. Detailed descriptions and results of 
the computation of the above cases are included as follows. 


A. 2-D Laminar Driven Square-Cavity Flows 

The first test case is concerning laminar recirculating flows inside a square 
cavity. Only one side of the walls is moving at a constant speed tangent to that 
wall. This case has been studied extensively by Burggraf [15] and has often been 
used as one of the standard testing cases for numerical methods in solving the 
incompressible Navier-Stokes equations. Physical geometry and wall boundary condi- 
tions are illustrated in Figure 4. Reynolds number of the flow (based on the cavity 
size and the moving wall velocity) studied in the present analysis is 400. Two 
different mesh systems, as shown in Figure 5, are used to study the effect of grid 
non-orthogonality on the accuracy of the present method. The grid system of Figure 
5(a) is uniform and orthogonal while the grid system of Figure 5(b) is non-uniform 
and non-orthogonal . 


UPPERWALL, v = o 
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Figure 4. Physical geometry and wall boundary conditions for 
laminar flows inside a wall-driven square cavity. 





Results of the computations are shown in Figures 6 and 7. Velocity vector 
plots of the predicted flow fields are compared in Figure 6 for the mesh systems 
shown in Figure 5. Detailed comparisons of the predicted velocity profiles along the 
mid- section of the cavity are illustrated in Figure 7. Predicted results of Burggraf 
[15] are also included. Good agreements between the present calculations and those 




Figure 6. Velocity vector plots, (a) Orthogonal grid, 
(b) Nonorthogonal grid. 



u ^upperwall 


Figure 7. Comparisons of velocity profiles along the mid- section 

of the square cavity. 
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of Burggraf [15] are also included. Good agreements between the present calcula- 
tions and those of Burggraf are shown in Figure 7. Discrepancies between the pre- 
sent predictions and Burggraf's results are mainly due to the hybrid differencing 
scheme used in the present method. The upwind part of the hybrid scheme produces 
large numerical diffusion which tends to reduce the strength of the vorticity inside 
the cavity. Effects of differencing schemes in approximating the convection terms on 
the predicted results will be studied in the next test case. 

Convergence history of the computation of the present case using uniform grids 
is given in Figure 8 which shows that the present numerical method is quite different. 
Almost identical convergence rates were found for the non-orthogonal case. 



_4 | I 1 1 

0 100 200 300 400 

ITERATION 


Figure 8. Convergence history for the driven cavity problem, Re = 400. 


B. 2-D Laminar Flows Over a Backward- facing Step 

This test case concerns 2-D laminar recirculating flows over a backward- facing 
step with 1:2 expansion ratio. The dependence of the size of the recirculation 
region (characterized by the reattachment length) on the Reynolds number (based 
on the inlet bulk velocity and twice of the inlet channel width) of the flow is of 
major concern. The physical domain and boundary conditions are illustrated in 
Figure 9 in which a fully developed laminar flow velocity profile is imposed at the 
flow entrance. A non-uniform grid of 45 x 45 was used for numerical computations. 
Several cases with different Reynolds numbers from 100 to 800 have been studied. 

An experimental and theoretical study about this problem, which results will be used 
as the basis of data comparisons, has been provided by Amaly et al. [16]. 
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dx 


o 
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dx 


o 


Figure 9. Physical geometry and boundary conditions of laminar flows 
over a backward-facing step (1:2 expansion). 

To save computational efforts, the solution of one case with Reynolds number 
100 is obtained in the first run. Then, a series of cases with increasing Reynolds 
numbers (i.e., 100, 200, 300, 400, 600, and 800) are calculated using the preceding 
results of lower Reynolds number as the initial guesses of the flow field. In this 
way , an average of 500 iterations for each case were needed to obtain converged 
solutions . 

Two different differencing schemes in approximating the convection terms are 
employed to demonstrate the effects of the differencing schemes on the predictions. 
One of the schemes is the widely used hybrid scheme [9], The other scheme employs 
the central differencing scheme plus an artificial dissipation term used to stabilize 
the solution which is similar to the one used by Rhie [17], The artificial dissipation 
term becomes effective only when the cell Peclet number (or cell Reynolds number) 
exceeds 10. 

Results of the present predictions using two different differencing schemes are 
compared with the experimental measurements [16] and other predictions as shown in 
Figure 10. It can be seen clearly from Figure 10 that the present method with hybrid 
scheme gives results similar to those predicted by TEACH code [16] while the present 
method with central differencing and artificial dissipation reveals predictions close to 
those predicted by INS 3D [18] and the method of Kim and Moin [19], This is reason- 
able since the TEACH code and the present method (with the first scheme) use the 
hybrid scheme which introduces large numerical dissipation by its upwind part (for 
cell Peclet number greater than 2) . This tends to reduce the reattachment length 
for Reynolds number greater than 400. The second scheme, which is similar to the 
ones used in INS 3D and the method of Kim and Moin, has the numerical accuracy 
close to second order by setting the artificial dissipation to be as small as the solu- 
tion stability permits such that better accuracy of the predictions is expected. 
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Figure 10. Reattachment length versus Reynolds number for laminar 
flows over a backward- facing step (1:2 expansion). 

Stream function plots of the predictions using the two differencing schemes for 
Reynolds number 600 are compared in Figure 11. It is shown in Figure 11 that the 
second scheme gives a smooth shape of the recirculation zone while the hybrid scheme 
gives a sudden change in the shape of the recirculation region upstream of the 
reattachment point. Also, larger sizes of the separation regions on the step side 
wall and along the upper wall are predicted using the second scheme. 




Figure 11. Streamline plots for laminar flow over a back ward- facing step 
(1:2 expansion), (a) Hybrid Scheme, (b) Central differencing 
plus artificial dissipation scheme. 
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C . 2-D Turbulent Flows Over a Backward-Facing Step 

In order to demonstrate the applicability of the present method to turbulent 
flow case, one of the standard test cases presented in the Stanford Conference [20] 
is selected here (i.e., turbulent flow over a 2:3 expansion backward- facing step). 
The standard k- e turbulence model was used to provide the eddy viscosity for the 
transport equations. The physical geometry and boundary conditions imposed are 
shown in Figure 12. The calculation domain extends upstream of the expansion plane 
by 4 step heights and downstream of the expansion plane by 30 step heights to 
assure a fully developed velocity profile at the exit. A uniform velocity profile is 
located at the inlet plane. A 45 x 42 grid was used in the computation. 300 itera- 
tions were required to obtain converged solutions. Only hybrid differencing schemes 
were used in this case. 


U = O, V = o 



du 

dx 

dv 

dx 


= o 


dx 


dk 

dx 

de 

dx 


= o 


= o 


Figure 12. Physical geometry and boundary conditions of turbulent flows 
over a backward- facing step (2:3 expansion). 

Results of the computation are shown in Figures 13, 14, and 15. These results 
are compared with the experimental measurements [20]. The under-prediction of the 
reattachment length is mainly due to the fast development of the mixing layer down- 
stream of the expansion plane which is the characteristics of the standard k-e tur- 
bulence model. Numerical diffusion provided by the hybrid scheme also contributes 
some part to the discrepancies between the predictions and measurements. 


D. Developing Laminar Flow Inside a 90-Deg-Bend Square Duct 

This test case simulates a three-dimensional developing laminar flow inside a 
90-deg-bend square duct as illustrated in Figure 16(a). The symmetry plane is 
located at z = 0 where the symmetric boundary conditions are imposed. A fully 
developed velocity profile of laminar flow inside a straight square duct is prescribed 
at the entrance which is 2.8 duct widths upstream of the bend, A zero pressure 
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PRESENT METHOD (k-e MOOEL) 

O EXP. [20] 



x/h 


Figure 13. Locus of flow reversal inside the recirculation region for 
turbulent flow over a backward- facing step (2:3 expansion). 



5.4h 

Figure 14. Stream line pattern of turbulent flow over a backward-facing 

step with 2:3 expansion ratio. 



2 

Figure 15. Contours of turbulent kinetic energy (k/U o ) of turbulent flow 
over a backward-facing step with 2:3 expansion ratio. 
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Figure 16. Geometry and mesh system of a 90-deg-bend square duct 

developing laminar flow problem. 

gradient exit (which is 4.5 duct widths downstream of the bend) boundary condition 
is imposed. The Reynolds number of the flow (based on the duct hydrolic diameter 
and the inlet bulk velocity) is 790. A 21 x 18 x 10 grid was used for numerical 
computations. The front view and side view of the mesh system are illustrated in 
Figure 16(b). Experimental measurements of Humphrey et al. [21] are used for data 
comparisons . 

Velocity vector plots on three sections along the main flow directions (i.e., on 
x-y plane) are shown in Figure 17. Secondary flow patterns at several stations 
across the bend are illustrated in Figure 18. These results are very similar to those 
obtained by Vanka [22] and Rhie [23]. Grid sizes of 50 x 22 x 15 and 58 x 15 x 11 
were used by Rhie and Vanka, respectively. The present investigation, using only 
less than half of their grid numbers, gives highly encouraging results. Detailed 
comparisons between the measured and the predicted main velocity profiles are given 
in Figure 19. 

With the above successful numerical simulations, it is believed that the present 
numerical method can be applied to general fluid dynamics problems with good 
numerical accuracy and efficiency . 
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Figure 17. Primary velocity patterns of laminar flow inside a 90-deg-bend square duct 

fal z/d = 0.0. (b) z/d = 0.25. (c) z/d = 0.48. 



Figure 18. Secondary velocity patterns of laminar flow inside a 90-deg bend, 
(a) 9 = 0 deg. (b) 6 = 30 deg. (c) 6 = 60 deg. (d) 9 = 90 deg. 






U/U r ef RHIE (22 x 15 x 50 GRID) [23] 

© EXP., HUMPHREY ET. AL. [21] 



u/u ref 

Figure 19. Primary velocity profiles for a 3-D 90-deg-bend 
square duct, (a) z/d = 0.0. (b) z/d = 0.25. 


CONCLUSIONS 


A numerical method for solving the steady or transient incompressible Navier- 
Stokes equations in three-dimensional body- fitted coordinate systems has been 
developed. In the present paper, the basic numerical algorithms and grid arrange- 
ments have been described in detail. A brief user’s guide to the present computer 
code (CNS3D) has been included in Appendix A. A program listing has also been 
attached in Appendix C . 

Several numerical testing examples of 2-D and 3-D, laminar and turbulent flow 
problems included in the present work have demonstrated that the present computer 
code is efficient and robust, and can be used as a reliable tool for engineering 
design and analysis applications. Applications of the present code to the internal 
turbulent flow problems of the SSME will be presented in the future publications. 
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APPENDIX A 


COMPUTER CODE STRUCTURES AND USER'S GUIDE 


The global structure of the present computer code (CNS3D) can be represented 
by a flow chart, shown in Figure A-l. The user is referred to Appendix C for 
detailed information. First, the program requests inputs, from logic unit 5 (LU = 5), 
of program control parameters that specify the maximum number of iterations, the 
type of flow (i.e., laminar or turbulent), number of iterations for solving the pres- 
sure correction equation (typically 10) , and underrelaxation factors for solving the 
transport equations, etc. This is followed by the definitions of all the program con- 
stants including turbulence model constants (these constants are subject to change 
according to the user's specific flow problem). Next, the program asks for inputs of 
the initial flow field guess from a restart file (LU = 8) which contains the grid 
system coordinates and flow field data that may be created by the user (including 
grid generation) or obtained from the previous solutions. Format of this data file is 
also subject to change according to the user's preference. Next, wall boundary con- 
trol parameters, boundary grid normal distance to the wall, and wall boundary direc- 
tion cosine are calculated in subroutine DIRCOS. Subroutine TRANF is then invoked 
to obtain the grid transformation coefficients. Before the solution procedure starts, 
the inlet mass flow rate is calculated which will be used to control the outlet mass 
flow rate to enhance mass conservation. The solution procedures consist of a series 
of subroutine calls to SOLVEQ starting from the solutions of the velocity vectors, u, 
v, and w, and then the solutions of scalar quantities (including the energy equation 
and the turbulence model equations) and finally the solution of the pressure correc- 
tion equation to update the velocity and pressure field such that a divergence- free 
flow field can be retained. 

After each global iteration of the solution procedures, the numerical of itera- 
tions and the maximum flow field corrections are checked with the initial settings. 

If the convergence criterion is satisfied or the number of iterations reaches the 
prescribed value then the solution procedures stop and the flow field solutions will 
be written on the pre-assigned disc file (LU = 7). 

For instance, if a steady- state laminar flow problem (Reynolds number of 600) 
is of interest and a converged solution is expected within 300 iterations and the 
number of iterations for solving the pressure correction equation is 10 and the 
underrelaxation factors are 0.5 and 0.95 for transport equations and pressure correc- 
tion equation, respectively, the first inputs from LU = 5 would be: 


Line 

1. 300 1 10 1 

2. 0.5 0.5 0.5 0.95 0.5 0.5 0.5 0.5 

3. 600. 0.0 

In the second input sequence (i.e., from restart file), the program reads in 
L x M x N lines of data records. See Figure A- 2 for grid structures. Notice that 
the program requires variable dimensions of (L+l, M+l, N+l) for solving the pressure 
correction equation. It is important to check the COMMON table for proper variable 
dimensions . 
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Figure A-l. Global structure of the present computer program CNS3D. 
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Figure A- 2. Grid mesh structures for 3-D calculations. 

If the flow problem involves symmetric or cyclic boundary conditions, then the 
user can look into the subroutine SYMOUT to specify the appropriate boundary con- 
ditions (the conditions shown in the program listing of Appendix C are for symmetric 
boundary conditions at K = 1). For cyclic boundary conditions at K = 1 and K = N, 
data at K = 2 and K = N-l can be used to obtain boundary conditions at K = 1 and 
K = N by requiring same gradients across K = 1 and K = N . This method is simple 
but will lag the boundary conditions by one iteration. A direct method without 
lagging the boundary conditions can also be employed by modifying the subroutine 
of linear algebra solver LINERX such that the boundary conditions can be part of the 
solution of the TDMA (tridiagonal matrix) solver. 

In case of incorporating different wall functions for turbulent flow problems 
(e.g.. References 24, 25, and 26), subroutine BOUNC and WALFN can be modified 
according to the user's method of wall treatments. The set of wall functions given 
in the program listing of Appendix C are derived from the conventional wall law and 
the equilibrium turbulent kinetic energy relations [8] . 

When additional source terms are to be added to the transport equations due to 
flow problem requirements, modifications to the source term calculation section in the 
subroutine SOLVEQ can be carried out. Notice that in the subroutine SOLVEQ source 
terms for the velocities v and w are included in the u-source section. Purpose of 
this is to save some computing time since these source terms use similar calculation 
routines . 

Some times it is required to solve more transport equations other than the basic 
ones included in Appendix C. To modify the program to incorporate more equations, 
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several changes are necessary. First, new variables must be added to the COMMON 
table (this can be easily done through the computer editor session). Then, new 
source term sections are added in the SOLVEQ subroutine. Finally, subroutines 
WALFN and SYMOUT are modified to incorporate the new variables into the boundary 
condition setting routines. 
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APPENDIX B 


A(K) 

AB(I, J,K) 

AE(I ,J,K) 

ALC 

ALE 

ALK 

ALP 

ALU 

ALV 

ALVIS 

ALW 

AN(I,J,K) 
ANAB 
ANVl(I) 
ANW1(I) 
AP(I , J ,K) 
APO(I ,J,K) 
ARDEN 
AREA 
AS(I ,J,K) 
AT(I ,J ,K) 
AW(I ,J,K) 
B(K) 

BB(I ,J,K) 
BOUNC 


LIST OF FORTRAN SYMBOLS 
= Matrix elements of a tridiagonal matrix 

= Link coefficients through the bottom face of a control volume 
= Link coefficients through the east face of a control volume 
= Underrelaxation factor for symmetry or cyclic boundary conditions 
= Underrelaxation factor for the e-equation 
= Underrelaxation factor for the k-equation 

= Underrelaxation factor for the pressure correction equation 
= Underrelaxation factor for the u-equation 
= Under relaxation factor for the v-equation 
= Underrelaxation factor for the effective viscosity 
= Underrelaxation factor for the w-equation 
= Link coefficient through the north face of a control volume 
= Sum of the link coefficients at all faces 
= Modified wall boundary link coefficient for v-equation 
= Modified wall boundary link coefficient for w-equation 
= Sum of the link coefficients around a control volume 
= Link coefficients in time marching direction 
= Area times density across a section in physical domain 
= Area of a section in physical domain 

= Link coefficients through the south face of a control volume 

= Link coefficients through the top face of a control volume 

= Link coefficients through the west face of a control volume 

= Matrix elements of a tridiagonal matrix 

= Coefficients in Stone’s partial factorization technique 

= Subroutine for getting turbulent wall boundary conditions through 
wall functions 



C(K) 

Cl 

C2 

CB 

CE 

CK 

CMU 

CMU1 
CMU 2 
CN 

cs 

CT 

CW 

CX(I , J,K) 
CY(I,J,K) 
CZ(I,J,K) 
D(K) 

DDB 

DDE 

DDN 

DDS 

DDT 

DDW 

DE(I , J ,K) 
DEO(I , J ,K) 
DEN (I , J ,K) 
DENO(I , J ,K) 
DENC 


= Matrix elements of a tridiagonal matrix 
= Turbulence model constant, = 1.44 
= Turbulence model constant, = 1.92 

= Convective flux through the bottom face of a control volume 

= Convective flux through the east face of a control volume 

= Von Karman constant, = 0.4 
= Turbulence model constant, = 0.09 

= CMU**0. 25 
= CMU**0. 75 

= Convective flux through the north face of a control volume 

= Convective flux through the south face of a control volume 

= Convective flux through the top face of a control volume 

= Convective flux through the west face of a control volume 

= Grid transformation coefficient, E 

s x 

= Grid transformation coefficient, E 

y 

= Grid transformation coefficient, E 

z 

= Matrix elements of a tridiagonal matrix 

= Diffusive flux through the bottom face of a control volume 

= Diffusive flux through the east face of a control volume 

= Diffusive flux through the north face of a control volume 

= Diffusive flux through the south face of a control volume 

= Diffusive flux through the top face of a control volume 

= Diffusive flux through the west face of a control volume 

= Turbulent kinetic energy dissipation rate, e 
= DE at the previous time level 
= Density of the fluid 
= DEN at the previous time level 
= Density at the center of a surface 
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DENIN 

DIRCOS 

DITM 

DK(I,J,K) 

DKO(I,J,K) 

DTT 

DU(I,J,K) 
DV (I , J ,K) 
DW(I,J,K) 

E 

EREXT 

ERRE 

ERRF 

ERRK 

ERRM 

ERRU 

ERRV 

ERRW 

EX(I,J,K) 

EY (I , J ,K) 

EZ(I,J,K) 

F(I,J,K) 

FO(I , J ,K) 

Fl(I.J.K) 

FLOW 

FLOWIN 

GEN (I , J ,K) 


= Initial value of density of the fluid 

= Subroutine for calculating the boundary grid sizes and direction 
cosines 

= Wall boundary average value of dissipation rate 
= Turbulent kinetic energy, k 
= DK at the previous time level 
= Time step size, At 

= Diffusive coefficient for the p’-equation 
= Diffusive coefficient for the p'-equation 
= Diffusive coefficient for the p’-equation 
= Wall law constant, = 9.01069 
= Convergence criterion tolerance 
= Maximum correction in e 
= Maximum correction of a variable 
= Maximum correction in k 
= Maximum correction in p 
= Maximum correction in u 
= Maximum correction in v 
= Maximum correction in w 
= Grid transformation coefficient, n x 
= Grid transformation coefficient, n 

y 

= Grid transformation coefficient, n 

z 

= Tentative variable of the transport equations 
= F at the previous time level 

= Variable quantity at the previous iteration step 
= Outlet mass flow rate 
= Inlet mass flow rate 

= Turbulent kinetic energy production rate 


33 



HINUM 

I 

IBC(I) 

IE 

IG 

IITO 

IITY 

IJLO(I,J,K) 

INIT 

INPRO 

INSOE 

INSOK 

INSOP 

INSOT 

INSOU 

INSOV 

INSOW 

IS 

ISWE 

ISWK 

ISWP 

ISWU 

ISWV 

ISWW 

IT 

ITT 

J 


= Large number, = 1.E30 
= Index along the £ grid lines 
= Boundary grid index 

= Index assigned for the transport equations 

= Problem control parameter, =1 for laminar flow and =2 for turbulent 
flow 

= Total number of wall boundary grids 
= Boundary grid face type 
= Boundary grid sequential order 
= Subroutine for initializing variables 
= Logical parameter for updating the effective viscosity 
= Logical parameter for solving the e-equation 
= Logical parameter for solving the k-equation 
= Logical parameter for solving the p '-equation 
= Logical parameter for solving the T- equation 
= Logical parameter for solving the u-equation 
= Logical parameter for solving the v- equation 
= Logical parameter for solving the w-equation 
= Starting value of I of the solution domain 
= Number of sweeps for solving the e-equation 

= Number of sweeps for solving the k-equation 

= Number of sweeps for solving the p'-equation 

= Number of sweeps for solving the u-equation 

= Number of sweeps for solving the v-equation 

= Number of sweeps for solving the w-equation 

= Last value of I of the solution domain 
= Number of time steps 
= Index along the n grid lines 


34 



JBC(I) 

JS 

JT 

K 

KBC(I) 

KS 

KT 

L 

LO 

LI 

L2 

LINERX 

LT 

M 

MO 

Ml 

M2 

MC(I,J,K) 

MT 

N 

NO 
N 1 
N2 

NEWVIS 

NLIMT 

NT 

P 

PCXI 


= Boundary grid index 

= Starting value of J of the solution domain 
= Last value of J of the solution domain 
= Index along the s grid lines 
= Boundary grid index 

= Starting value of K of the solution domain 
= Last value of K of the solution domain 
= Maximum dimension of grid system in I direction 
= L + 1 

= Starting point of blockage region in I direction 
= Last point of blockage region in I direction 
= Subroutine for solving algebraic equations 
= L - 1 

= Maximum dimension of grid system in J direction 
= M + 1 

= Starting point of blockage region in J direction 
= Last point of blockage region in J direction 
= Wall blockage region control parameter 
= M - 1 

= Maximum dimension of grid system in K direction 
= N + 1 

= Starting point of blockage region in K direction 
= Last point of blockage region in K direction 
= Subroutine for updating the effective viscosity 
= Limit of maximum number of iterations 
= N - 1 

= Static pressure (relative) 

= Pressure gradient, P^ 



PDUV 

PEDA 

PP 

PPBLK 

PSCI 

PTA 

PW 

RENL 

SIGE 

SIGK 

SINX(I) 

SINY(I) 

SINZ(I) 

SMNUM 

SOC1 

SOC2 

SOC3 

SOLVEQ 

SP(I,J,K) 

SPK(I,J,K) 

SU(I,J,K) 

SUK(I , J ,K) 

SX(I , J ,K) 

SY(I,J,K) 

SYMOUT 

SZ(I,J,K) 

TAUN(I) 


= Blockage control parameter for link coefficients 
= Pressure gradient, 

= Pressure correction, p* 

= Global pressure correction 
= Pressure gradient, P^ 

= Wall boundary source term for the momentum equations 
= Wall value control parameter 
= Reynolds number of the fluid 
= Turbulence model constant, = 1.3 
= Turbulence model constant, = 1.0 
= Wall boundary direction cosine 
= Wall boundary direction cosine 
= Wall boundary direction cosine 
= Small number, l.E-30 
= Source term due to shear stress 

= Source term due to shear stress 

= Source term due to shear stress 

= Subroutine for solving general transport equation 
= Linear part of the source term 
= Secondary linear part of the source term 
= Constant part of the source term 
= Secondary constant part of the source term 
= Grid transformation coefficient, c 
= Grid transformation coefficient, £y 
= Subroutine for setting flow boundary conditions 
= Grid transformation coefficient, r, 

z 

= Wall shear stress 
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TIMT 

TJO(I,J,K) 
TM(I,J,K) 
TMO(I ,J,K) 
TMULT 
TRANF 
TXXE(I,J,K) 
TXXW(I,J,K) 
TXYN (I , J ,K) 
TXYS(I,J,K) 
TXZT(I,J,K) 
TXZB(I , J ,K) 
TYYN (I , J ,K) 
TYYS(I , J ,K) 
TYXE(I ,J,K) 
TYXW(I , J , K) 
TYZT(I ,J,K) 
TYZB (I , J ,K) 
TZZT (I , J ,K ) 
TZZB (I , J ,K) 
TZXE(I , J ,K) 
TZXW(I,J,K) 
TZYN (I , J ,K) 
TZYS(I,J,K) 
U(I,J,K) 
UO(I , J,K) 

UC 

UCXI 


= Total time 

= Jacobian of metric transformation 
= Temperature 

= TM at the previous time level 
= Wall shear stress 

= Subroutine for calculating the grid transformation coefficients 
= Metric coefficient for east face diffusive flux 
= Metric coefficient for west face diffusive flux 
= Metric coefficient for north face diffusive flux 
= Metric coefficient for south face diffusive flux 
= Metric coefficient for top face diffusive flux 
= Metric coefficient for bottom face diffusive flux 
= Metric coefficient for north face diffusive flux 
= Metric coefficient for south face diffusive flux 
= Metric coefficient for east face diffusive flux 
= Metric coefficient for west face diffusive flux 
= Metric coefficient for top face diffusive flux 
= Metric coefficient for bottom face diffusive flux 
= Metric coefficient for top face diffusive flux 
= Metric coefficient for bottom face diffusive flux 
= Metric coefficient for east face diffusive flux 
= Metric coefficient for west face diffusive flux 
= Metric coefficient for north face diffusive flux 
= Metric coefficient for south face diffusive flux 
= U- velocity 

= U at the previous time level 
= Velocity at the center of a surface 
= U- velocity gradient, u^ 



UEDA 

= U- velocity gradient, 

UINC 

= Velocity correction at outlet plane 

USCI 

= U- velocity gradient, 

UX 

= U- velocity gradient, u x 

UY 

= U- velocity gradient, Uy 

UZ 

= U- velocity gradient, u 

z 

V (I , J , K) 

= V- velocity 

VO(I,J,K) 

= V at the previous time level 

vise 

= Molecular viscosity, y 

VISE(I,J,K) 

= Effective viscosity, y^ 

VCXI 

= V- velocity gradient, v 

VEDA 

= V -velocity gradient, v^ 

VSCI 

= V- velocity gradient, v^ 

VX 

= V-velocity gradient, v x 

VY 

= V-velocity gradient, v 

VZ 

= V-velocity gradient, v 

z 

W(I , J ,K) 

= W- velocity 

WO(I,J,K) 

= W at the previous time level 

WALLFN 

= Subroutine for calculating the wall functions 

WALVAL 

= Subroutine for assigning wall values 

WCXI 

= W- velocity gradient, w^ 

WEDA 

= W- velocity gradient, w^ 

WSCI 

= W- velocity gradient, w^ 

WX 

= W- velocity gradient, w x 

WY 

= W- velocity gradient, w 

WZ 

= W-velocity gradient, w 

z 

X(I , J ,K) 

= X -coordinate 
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Y(I,J,K) 

= Y -coordinate 

YN(I) 

= Wall normal distance from the last grid 

YN 1(1) 

= Wall grid volume size 

YPLN(I) 

= Non dimen sionalized YN, y + = u^y/v 

Z(I,J,K) 

= Z -coordinate 
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